/*
 * SpeciesUtils.cpp
 *
 *  Created on: 22 Aug 2011
 *      Author: Allan
 */

#include "SpeciesUtils.h"

// C++ includes
#include <set>

// Boost includes
#include <boost/foreach.hpp>

// GeoReact includes
#include "Geochemistry/GeoDatabase.h"

#include "Utils/Algorithms.h"
#include "Utils/StringUtils.h"

const vector<string> FilterIons(const vector<string>& species)
{
	// The list of aqueous ions from the given list of species
	vector<string> ions;
	
	// Loop over all the species and check which of them is an aqueous ion
	BOOST_FOREACH(const string& s, species)
		if(s.find("+") != s.npos || s.find("-") != s.npos)
			ions.push_back(s);
	
	return ions;
}

const vector<string> FilterComplexes(const vector<string>& species)
{
	// Get an instance of the database
	const GeoDatabase& db = GeoDatabase::GetInstance();
	
	// The list of aqueous complexes from the given list of (species)
	vector<string> complexes;
	
	// Loop over all the species and check which of them is an aqueous complex
	BOOST_FOREACH(auto s, species)
		if(db.FindComplexInfo(s))
			complexes.push_back(s);
	
	return complexes;
}

const vector<string> FilterBasisSpecies(const vector<string>& species)
{
	return FilterSpecies(species, {"basis"});
}

const vector<string> FilterSecondarySpecies(const vector<string>& species)		
{
	return RemoveSpecies(species, {"basis"});
}

const vector<string> FilterSpecies(const vector<string>& species, const vector<string>& tags)
{
	// The subset of the given species that have any of the given tags
	vector<string> newSpecies;
	
	// Loop over all the given species 
	BOOST_FOREACH(const string& s, species)
		if(GeoDatabase::GetInstance().GetSpeciesInfo(s).HasAnyTag(tags))
			newSpecies.push_back(s);
			
	return newSpecies;
}

const vector<string> RemoveSpecies(const vector<string>& species, const vector<string>& tags)
{
	// The subset of the given species that does not have any of the given tags
	vector<string> newSpecies;
	
	// Loop over all the given species 
	BOOST_FOREACH(const string& s, species)
		if(!GeoDatabase::GetInstance().GetSpeciesInfo(s).HasAnyTag(tags))
			newSpecies.push_back(s);
			
	return newSpecies;
}

const double ElectricalCharge(const string& species)
{
	vector<string> words = Split(species, "[]");
	
	if(words.size() == 1) return 0;
	
	if(words[1] == "+") return +1;
	if(words[1] == "-") return -1;
	
	const string strCharge = {words[1][1], words[1][0]};
	
	const double charge = atof(strCharge.c_str());
	
	return charge;
}

const double ElementAtomsNumeral(const string& species, unsigned index)
{
	// This ensures the distinction between the elements Fe and F, Ni and N, Sn and S, etc.
	if(islower(species[index]))
		return 0;
	
	if(index >= species.size())
		return 1;
	
	unsigned size = 0;
	
	for(; isdigit(species[index + size]); ++size);
	
	if(size == 0) 
		return 1;
	else
		return atof(species.substr(index, size).c_str());
}

const double ElementAtomsFactor(const string& species, unsigned index)
{
	const unsigned i1 = species.find(')', index);
	const unsigned i2 = species.find('(', index);
	
	if(i1 < i2)
		return ElementAtomsNumeral(species, i1+1);
	else
		return 1;
}

const double ElementAtomsFormula(const string& species, const string& element)
{
	unsigned pos = 0;
	
	double n = 0;
	
	while(true)
	{
		pos = species.find(element, pos);
		
		if(pos == unsigned(-1))
			break;
		
		pos += element.size();
		
		n += ElementAtomsNumeral(species, pos) * ElementAtomsFactor(species, pos);
	}
	
	return n;
}

const double ElementAtoms(const string& species, const string& element)
{
	optional<ReactionInfo> info = GeoDatabase::GetInstance().FindReactionInfo(species);
	
	if(!info) return ElementAtomsFormula(species, element);
	
	double n = 0.0;
	
	for(unsigned i = 0; i < info->products.size(); ++i)
	{
		n += info->stoichiometries[i] * ElementAtomsFormula(info->products[i], element);
	}
	
	return n;
}
